// eee.Sheffield.PZ.Imaging
//
// Copyright ?Ping Zou, 2007
// sg71.cherub@gmail.com

using System;
using System.Collections.Generic;
using System.Text;
using System.IO;
using eee.Sheffield.PZ.Math;

namespace eee.Sheffield.PZ.Imaging.ImageFilter
{
    public enum ProfileShape { Valley, Peak };

    // Sun 1989
    public class SunMatchedFilter
    {
        public PZMath_vector RectangularMatchedFilter(PZMath_vector sl, int initr)
        {
            //int width = sl.Size;
            //int outputWidth = width - initr > 0 ? width - initr : 0;
            //PZMath_vector output = new PZMath_vector(outputWidth);

            //int lowerBound = initr < outputWidth ? initr : 0;
            //int upperBound = outputWidth - 1;

            //if (lowerBound >= upperBound)
            //    throw new ApplicationException("lowerBound should be less than upperBound");

            //for (int i = lowerBound; i <= upperBound; i++)
            //{
            //    double iterm1 = 0.0;
            //    double iterm2 = 0.0;
            //    double iterm3 = 0.0;
            //    for (int j = (i - initr); j <= i + initr; j++)
            //        iterm1 += sl[j];
            //    for (int j = 0; j <= i - initr - 1; j++)
            //        iterm2 += sl[j];
            //    for (int j = (i + initr + 1); j <= width - 1; j++)
            //        iterm3 += sl[j];
            //    output[i] = iterm1 - iterm2 - iterm3;
            //}

            int width = sl.Size;
            // generate rectangular matched filter
            PZMath_vector m = new PZMath_vector(width);
            for (int i = 0; i < width; i++)
            {
                if (i < initr)
                    m[i] = -1.0;
                else if (i < 2.0 * initr)
                    m[i] = 1.0;
                else
                    m[i] = -1.0;
            }

            // do convolution
            PZMath_vector output = m.Convolve(sl);

            //output.WriteFile("con_1.txt");
            //int max = output.MaxIndex();

            return output;
        } // RectangularMatchedFilter()

        public PZMath_vector TriangularMatchedFilter(PZMath_vector sl, int initr)
        {
            int width = sl.Size;
            double halfWidth = (double) width / 2;
            // generate triangular matched filter
            PZMath_vector m = new PZMath_vector(width);
            for (int i = 0; i < width; i++)
                m[i] = 1.0 - 2 * System.Math.Abs(i - halfWidth) / (double)width;
            
            // do convolution
            PZMath_vector output = m.Convolve(sl);

            return output;
        } // TriangularMatchedFilter()

        public double SignalLevel(PZMath_vector sl, double c, double initr)
        {
            int width = sl.Size;
            // get signal region [c - r, c + r];
            int lowerBound = c - initr > 0 ? (int)(c - initr) : 0;
            int upperBound = c + initr < width ? (int)(c + initr) : width - 1;

            // calculate signal level
            double S = 0.0;
            for (int i = lowerBound; i <= upperBound; i++)
                S += sl[i];
            return S / (double)(upperBound - lowerBound + 1);
        } // SignalLevel()

        public double BackgroundLevel(PZMath_vector sl, double c, double initr)
        {
            int width = sl.Size;
            // get background region [0, c - r - 1] && [c + r + 1, width - 1]
            int lowerBound1 = 0;
            int upperBound1 = c - initr - 1 > 0 ? (int)(c - initr - 1) : 0;
            upperBound1 = upperBound1 < width ? upperBound1 : width;

            int lowerBound2 = c + initr + 1 < width ? (int)(c + initr + 1) : width - 1;
            lowerBound2 = lowerBound2 > 0 ? lowerBound2 : 0;
            int upperBound2 = width - 1;
            //int length = upperBound2 - lowerBound2 + 1 + upperBound1 - lowerBound1 + 1;

            // calculate background level
            int count = 0;
            double B = 0.0;
            for (int i = lowerBound1; i <= upperBound1; i++)
            {
                B += sl[i];
                count++;
            }
            for (int i = lowerBound2; i <= upperBound2; i++)
            {
                B += sl[i];
                count++;
            }

            return B / (double)count;
        } // BackgroundLevel()

        public double RollOffPoint(double S, double B)
        {
            return (S + B) / 2;
        } // RollOffPoint()

        public PZMath_vector EdgePoints(PZMath_vector sl, double rollOffPoint, ProfileShape profileShape)
        {
            PZMath_vector edgePoints = new PZMath_vector(2);
            
            int width = sl.Size;
            // find the left edge
            for (int i = 0; i < width; i++)
            {
                if (profileShape == ProfileShape.Peak)
                {
                    if (sl[i] > rollOffPoint)
                    {
                        edgePoints[0] = (double)i;
                        break;
                    }
                }
                else
                {
                    if (sl[i] < rollOffPoint)
                    {
                        edgePoints[0] = (double)i;
                        break;
                    }
                }
            }
            // find the right edge
            for (int i = width - 1; i > edgePoints[1]; i--)
            {
                if (profileShape == ProfileShape.Peak)
                {
                    if (sl[i] > rollOffPoint)
                    {
                        edgePoints[1] = (double)i;
                        break;
                    }
                }
                else
                {
                    if (sl[i] < rollOffPoint)
                    {
                        edgePoints[1] = (double)i;
                        break;
                    }
                }
            }

            return edgePoints;
        } // EdgePoints()  

        public void Sun(double initr, PZMath_vector sl,
            out double radius, out double centre, out double leftEdge, out double rightEdge)
        {
            // Sun's matched filter
            double sun_estimateR = initr;
            PZMath_vector sun_output = RectangularMatchedFilter(sl, (int)sun_estimateR);
            // find the maximum output position
            int sun_m = sun_output.MaxIndex() - sl.Size / 2;
            // get signal level and backgroud level
            double sun_signalLevel = SignalLevel(sl, sun_m, sun_estimateR);
            double sun_backgroundLevel = BackgroundLevel(sl, sun_m, sun_estimateR);
            // get roll-off point
            double sun_rollOffPoint = RollOffPoint(sun_signalLevel, sun_backgroundLevel);
            // find edge points
            PZMath_vector sun_edgePoints = EdgePoints(sl, sun_rollOffPoint, ProfileShape.Peak);
            // update radius and centre point 
            double sun_r = (sun_edgePoints[1] - sun_edgePoints[0] + 1.0) / 2.0;
            double sun_centre = sun_edgePoints[0] + sun_r;

            // output
            radius = sun_r;
            centre = sun_centre;
            leftEdge = sun_edgePoints[0];
            rightEdge = sun_edgePoints[1];
        } // Sun()

        public void SunRecord(double initr, PZMath_vector sl, string fileName)
        {
            FileStream stream = new FileStream(fileName, FileMode.Create, FileAccess.Write);
            StreamWriter writer = new StreamWriter(stream);

            // writer sl
            int slLength = sl.Size;
            for (int i = 0; i < slLength; i++)
            {
                writer.Write(String.Format("{0, -10:f}", sl[i]));
            }
            writer.WriteLine();

            // Sun's matched filter
            double sun_estimateR = initr;
            PZMath_vector sun_output = RectangularMatchedFilter(sl, (int)sun_estimateR);

            // write matched filter output
            int outputLength = sun_output.Size;
            for (int i = 0; i < outputLength; i++)
            {
                writer.Write(String.Format("{0, -10:f}", sun_output[i]));
            }
            writer.WriteLine();


            // find the maximum output position
            int sun_m = sun_output.MaxIndex() - sl.Size / 2;
            // writer maximum output position
            writer.WriteLine(sun_m);

            // get signal level and backgroud level
            double sun_signalLevel = SignalLevel(sl, sun_m, sun_estimateR);
            double sun_backgroundLevel = BackgroundLevel(sl, sun_m, sun_estimateR);
            // get roll-off point
            double sun_rollOffPoint = RollOffPoint(sun_signalLevel, sun_backgroundLevel);
            // find edge points
            PZMath_vector sun_edgePoints = EdgePoints(sl, sun_rollOffPoint, ProfileShape.Peak);
            // writer edge points
            writer.WriteLine(String.Format("{0, -10:f}{1, -10:f}", sun_edgePoints[0], sun_edgePoints[1]));

            // update radius and centre point 
            double sun_r = (sun_edgePoints[1] - sun_edgePoints[0] + 1.0) / 2.0;
            double sun_centre = sun_edgePoints[0] + sun_r;
            writer.WriteLine(String.Format("{0, -10:f}", sun_centre));

            writer.Flush();
            writer.Close();
            stream.Close();
        } // SunRecord()
    }
}
